Lattice-switch Monte Carlo for binary hard-sphere crystals 
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We show how to generalize the Lattice Switch Monte Carlo method to calculate the phase diagram 
of a binary system. A global coordinate transformation is combined with a modification of particle 
diameters, enabling the multi-component system in question to be explored and directly compared 
to a suitable reference state in a single Monte Carlo simulation. We use the method to evaluate 
the free energies of binary hard sphere crystals. Calculations at moderate size ratios, a = 0.58 and 
a = 0.73, are in agreement with previous results, and confirm AB2 and AB13 as stable structures. 

£ — We also find that the AB(CsCl) structure is not entropically stable at the size ratio and volume at 

which it has been reported experimentally, and therefore that those observations cannot be explained 

f^**) by packing effects alone. 

<N 
W) 

I. INTRODUCTION 

CO 

The most efficient packing of hard spheres is one of the most easily posed questions in physics. For monodisperse 
spheres, it forms the basis of Kepler's conjecture[TJ H], and the solution is any one of many degenerate close-packing 
arrangements. In physical systems such as opals and colloidal suspensions, the problem is generalised to calculating the 
highest entropy structure of a given density. Away from the degenerate close-packing limit, face-centred cubic turns 
out to be the most stable, and its exact entropy advantage has been accurately calculated [3J. A further generalization, 
which we consider here, is the packing of spheres of different radii: the so-called binary hard sphere system. 

In order to calculate stable packings for binary hard spheres, one must evaluate the free energies of candidate 
phases as a function of density, composition and size ratio, so that the entire phase diagram can be constructed by 
| ', minimization of the total free energy. 

A common approach to free energy calculation is use thermodynamic integration. This involves the construction 
of a reversible path between the structure of interest and some reference system. By performing a simulation or set 
of simulations that traverse that path, one can estimate the free energy difference between the reference and target 
systems by integrating the derivative of the free energy along that path. For example, the integration technique of 
Frcnkcl and Ladd [4 uses an Einstein crystal as a reference system, and has been successfully applied to a wide 
range of target systems, including binary hard sphere solids [3J EJ EJ |7] . However, the integration process introduces 
an intrinsic source of systematic error. Small systematic errors can accumulate at each step in the integration, and 
integrating across the singularity at a phase transition causes problems. Of course, careful application can minimise 
these effects, but it remains preferable to have a technique where systematic errors are more easily quantifiable. The 
most accurate free energy differences in hard sphere systems are obtained from lattice-switch Monte Carlo, [3] which 
£^ allows two phases to be sampled directly during a single Monte Carlo simulation, mapping between them using a global 
■^J- coordinate transformation. Provided some biased sampling method ensures that the phase transformation move is 
accepted, the free-energy difference between phases can be determined directly from the equilibrium probability of 
qq finding the simulation in each phase. Finite size effects are the only source of systematic errors. This approach was 
first used to study the solid state phase behavior of monodisperse hard spheres [51 H]. The method was subsequently 
t — ■ developed for different energy models: "Hamiltonian Switch" and off-lattice: "phase-switch". It has been applied to 
the freezing transition for hard spheres in [3|TU] and also to soft potentials in (TTJ [T5J [T3J . 

Here we investigate whether this methodology can be applied successfully to multi-component systems by general- 
ising the lattice-switch approach to examine the structural phase behaviour of binary hard spheres. The binary hard 
sphere system is well studied, and therefore the validity of the generalised lattice switch method can be verified by 
direct comparison with the results from previous experimental and theoretical work. 

Stable phases of the binary-hard sphere system are characterised primarily by the diameter ratio, a — obI^a-, 
where the large and small particles are designated as A and B respectively. At high values of a, from 1.0 down to 
« 0.92, the system is expected to form a substitutionally disordered crystal, where A and B particles are randomly 
distributed on an FCC lattice [2J [TS]. At slightly lower values, in the range 0.92 > a > 0.85, this mixed crystal 
becomes metastable, and the solid phase phase separates into FCC crystals of A and B particles pT4l IT5] . 
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Figure 1: The three binary structures considered in this work, with the A particles shown as pale spheres, B particles shown 
as darker spheres, and the simulation cell shown as a translucent box. a) The AB13 structure consists of a simple cubic lattice 
of large spheres with icosahedral clusters of 13 B particles in the body-center of each cell. The 112 particle unit cell contains 
8 B-clusters arranged in alternating orientations, where the clusters are rotate by 90° between adjacent sub-cells, b) In AB2, 
hexagonal close-packed layers of A particles are stacked directly on top of each other, with a honeycomb pattern of B particles 
places in the interstitial sites between the A layers, c) The AB(CsCl) structure, consisting of two interlaced simple cubic 
lattices, with the B particles at the body-centre of the A cube. 

Below this, for 0.85 > a > 0.5, three different lattices (shown in Figure[l} have been observed experimentally: AB2, 
AB13 and AB(CsCl). AB13 and AB2 have been observed in Brazilian gem opals (colloidal silica crystals) [T6| I17j. 
and later in poly-methylmethacylate (PMMA) colloidal dispersions, both of which closely resemble hard spheres. In 
PMMA, AB13 was first seen in [IS] for a « 0.61, followed by AB2 and AB13 at a m 0.58 in [TS]. These two structures 
were determined to be unstable above a w 0.62 in [20] and below a w 0.52 in [21]. More recent experiments have 
found AB2 to be stable for 0.60 > a > 0.425 and AB13 stable over 0.62 > a > 0.485 [22]. 

These complex lattices have also seen in systems that bear little resemblance to hard-spheres, such as rare gas solids 
[23], nanocrystals [23], charge-stabilised colloids [25], and block copolymer micelles [23]. At first, it is surprising that 
these structures, particularly the highly complex AB13 supcrlattice, should form spontaneously in such a wide range 
of systems. However, it appears that entropic effects of hard-sphere packing are sufficient to stabilise these structures. 

Previously, the lattice switch approach has been used for the free energy difference between two candidate structure 
for a single species. Two-species phase behaviour is more complicated, as the stable state of a binary system may 
consist of co-existing phases of different compositions, at equal pressures and chemical potentials. Therefore, instead 
of just comparing free energy of pairs of candidate phases, we must to evaluate the absolute free energies of all 
competing phases so that the overall phase diagram can be constructed. The basic approach has been presented in 
[27] by Bartlett, who considered only the fluid and pure (face-centered cubic A or B) crystal phases. We shall follow 
the work of Eldridge et al [5] [6] [2D] and Cottin and Monson (28] [29] and also include AB2, AB13 and AB(CsCl) as 
candidate two-component crystals. In the Eldridge [5] work, the stability of AB2 and AB13 for binary hard spheres 
was determined using the thermodynamic integration method of Frenkel and Ladd [3], with further details of the 
AB2 calculations in [6] and AB13 in [7j. These results show general agreement with the experimental systems, but 
suggest that AB2 is kinetically repressed at the AB2 1:2 stoichiometry, as it only appears at higher concentrations of 
B particles [22] . 

Theoretical studies of AB2, AB13 & AB(CsCl) have been also carried out via cell theory [25] [23], and density 
functional theory (DFT) [30, 31, 32 . For AB2 and AB13, the results arc generally consistent with both the integration 
method and experimental results, although some minor differences remain. In particular, the range of stability of 
AB13 is computed via cell theory to be 0.61 > a > 0.54 [23], whereas the integration method finds AB13 to be stable 
over the range 0.626 > a > 0.474 [7J. Although the latter work does not consider the possibility that AB2 is more 
stable than AB13 over this range of a, it is notable that the range of stability determined by integration methods is 
consistent with the range of observation of AB13 found in the experimental work presented in [22] , 

The case of the AB(CsCl) structure is somewhat less clear. The close-packed density of this structure peaks at 
a = 0.732, and so may compete with the slightly denser FCC structure, and indeed metastable CsCl has been observed 
experimentally of at a — 0.736 in [33] , The stability of this structure has not been ascertained by integration methods: 
it was not considered to be a valid candidate structure in [20 . The cell theory work in [25] did not find CsCl to be 
stable for any a, but DFT results have shown metastablity [31], and even stablity [32]. However, the work in [32] did 
not consider full phase separation into coexisting A and B FCC crystals; they only consider a fluid co-existing with a 
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crystal phase made up of randomly distributed A and B particles, but rich in A, and the fluid composition appears to 
be fixed at 1:1. Interestingly the AB(CsCl) structure has been stabilised experimentally by using oppositely charged 
colloids [25J|2H|351[3B]. Moreover, BCC is the lowest energy arrangement for charged monodisperse particles, so even 
if the AB(CsCl) structure is only metastable in general, it is possible that a small amount of charge transfer between 
particles or between particles and fluid is enough to stabilise the structure. 

For even smaller sphere size ratios (a < 0.44), the situation becomes more complicated as additional structures start 
to appear (e.g. AB(rocksalt) [37]) and because, at lower volume fractions (or sufficiently low a), the small spheres are 
free to move from site to site within the crystal of larger spheres. Such systems are not well suited to examination via 
the lattice-switch approach, as the algorithm would have to be modified heavily in order to cope in the free movement 
of small spheres. At smaller size ratios still, the success of the perturbative approach [38] suggests that it is best to 
model the presence of very small spheres via a depletion potential. 

Given the range of results that exist for binary hard spheres, it provides an excellent testing ground for the extension 
of the lattice-switch approach to multi-component systems. As our primary concern is simply to extend the lattice- 
switch technique to the binary system, it makes sense to choose a particular value of a to examine, instead of trying 
to explore a whole range of diameter ratios before the correctness of the approach has been determined. We have 
chosen to examine the stability of AB2 and AB13 at a = 0.58, as there has been a large amount of theoretical and 
experimental work published for that particular diameter ratio. Furthermore, due to the ambiguity concerning the 
stability of the AB(CsCl) structure, we shall also use the lattice-switch approach to examine the stability of AB(CsCl) 
relative to the fluid phase and phase-seperated solid state at a = 0.73. 



II. FORMULATION 



We consider a system of N particles, of spatial coordinates {f} and diameters {c}, confined within a volume V, and 
subject to periodic boundary conditions. The configurational energy for this system of hard spheres has the form 



E({r}) = 



if rij < ay V 
oo otherwise, 



(1) 



where = \fl — fj\ and = | (er, + (Xj). Each microstate has a Boltzmann weight exp[— 0E], and so the total 
configurational weight associated with this system can be written as 



m, to = n 17 H n e ^ - ^ 



(2) 



where Q(x) = 1(0) for x > 0(< 0), and the product on (ij) extends over all particle pairs. The associated entropy 
density is 



s(N,V) = -lne(N,V). (3) 

This formulation describes the total entropy, integrating over all possible positions and so over all possible phases. 
In order to compare the statistical weights of the configurations associated with individual (candidate) phases, we 
must formulate a constraint that identifies a configuration as 'belonging to' a given phase. To do this, we decompose 
the particle position coordinates into a sum of 'lattice' and 'displacement' vectors: 



f t = Rf + u l: (4) 

where {R}a = — 1, JV" forms a set of fixed vectors associated with the crystal structure labeled A. This 

set is defined as the orthodox crystallographic lattice convolved with the orthodox basis, but will be referred to here 
as simply the 'lattice vectors'. In previous work the coordinates have denoted a pure phase arrangement of the 
particles. However, it is possible for A to represent a two-phase mixture or for particle sizes to be changed in the 
switch. We are thus able to initialise a simulation within a particular phase, using the appropriate set of lattice vectors 
({R}a) an d diameters ({c}^) for that phase, just by setting iZj = V i. We can now use a simple spatial criterion 
to decide whether a particular configuration belongs to the same phase as the simulation progresses. Periodically, 
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throughout the simulation, the particles are checked to determine whether, after having taken any centre-of-mass drift 
into account, Ui is comparable to the interparticle spacing (for more details, see section III C I . This can only happen 



if the particles have been rearranged and the structure has been modified, so provides a robust criterion for ensuring 
that a set of microstates belong to the same macrostate. We note that this restriction of the configurational integral 
to pure perfect phases erases contributions from equilibrium concentrations of lattice defects, for example lattice 
vacancies, however for the systems in question such defects are extremely rare and make no significant contribution 
to total free energy [35] . 

By performing this coordinate decomposition (eq. [4]), we have effectively changed our stochastic variables from 
being {r} to [{it}, .A], operating under a modified configuration energy function (c.f. eq. [T| of the form 



E({u},A) = i ""'I;' 

' 1 oo otherwise, 



(5) 



where the lattice label A dictates the choice of lattice vectors via ry {{u}, A) = (Rf + Ui) — (Rf + Uj) 



diameters via <Jij{A) 
written as 



and 



aij) . The configuration weight associated with a particular structure can now be 



(6) 



where f. signifies integration subject to the single-phase constraint. The associated entropy density becomes 



s(N,V,A)=-lnn(N,V,A) 



(7) 



Once we construct a simulation capable of visiting the two different phases {A, B) , the entropy difference between 
them can be determined as: 



As AB = s(N, V, A) ~ s(N, V, B) = i \nTl AB {N, V), (8) 



where 



n(N,V,A) _ P(A\N,V) 

Uab{n > v) ~ q(n,v,b) - p(b\n,vy (9) 

Here P(A\N, V) is the joint probability that a system free to visit both structures will be found in the A structure. 



III. IMPLEMENTATION AND METHODOLOGY 



A. Defining the structures and the switch 



To evaluate the free energy of a particular binary crystal, we need to perform a lattice switch simulation capable of 
visiting the target phase, A, and some reference phase, B. As an accurate semi-empirical expression for the free energy 
of a monodisperse hard sphere FCC crystal has been determined by Alder et al SI] , we choose this structure to be 
our usual reference state. A lattice switch which kept the diameters of the particles fixed (such that {cr}.4 = 
would require either a single simulation box containing co-existing crystals of A and B particles (and thus associated 
interfacial effects, making for a poor reference system) or two separate boxes corresponding to separate A and B 
crystals. This latter approach can be made to work, but the complexity associated with adding a second simulation 
cell is unnecessary when the two cells are being used to simulate two identical structures that differ only in scale (by 
a). A simpler option is to allow the diameters to change during the switch, turning all o~b particles in the mixed 
phase into particles of diameter a a, allocated different sites on a single shared FCC lattice. This creates a switch 
between a binary crystal and a single crystal, and ensures that the number of degrees of freedom in each system is the 
same. Furthermore, as the free-energy difference between FCC and HCP is accurately known [S], we are free to use 
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Target Crystal 


Reference Crystal 


Runtime 


Runtime 


Runtime 


N 


Structure 


Dimensions 


Structure 


Dimensions 


MCS [EOS] 


MCS [WF] 


MCS [DF] 


112 


AB13 


lxlxl 


HCP 


4x7x4 


5 x 10 5 


10 6 [40] 


4 x 10 6 [8] 


896 


AB13 


2x2x2 


FCC 


8 x 8 x 14 


5 x 10 5 


6 x 10 6 [240] 


6 x 10 6 [10] 


192 


AB2 


4x4x4 


FCC 


8x4x6 


10 6 


10 6 [40] 


4 x 10 6 [80] 


648 


AB2 


6x6x6 


FCC 


8x9x9 


10 6 


5 x 10 6 [200] 


2 x 10 6 [200] 


54 


CsCl 


3x3x3 


FCC 


6x3x3 


5 x 10 5 


10 6 [40] 


4 x 10 6 [40] 


128 


CsCl 


4x4x4 


HCP 


8x4x4 


5 x 10 5 


2 x 10 6 [80] 


5 x 10 6 [20] 


432 


CsCl 


6x6x6 


FCC 


8x6x9 


5 x 10 5 


3 x 10 6 [120] 


8 x 10 6 [40] 



Table I: Structures and sizes used in this work. For FCC, HCP and AB2, the unit-cells were oriented so that the close-packed 
layers were lying in the x-y plane. The cubic cells of AB13 and CsCl were made commensurate with the simulation cell. The 
run times are shown for determining the equation of state (EOS), the weight function (WF) and the free-energy difference 
(DF), in units of Monte Carlo sweeps (MCS). Numbers in square brackets indicates the largest number of blocks the run was 
split into for block analysis. For all structures and sizes, 20,000 MCS runs were used to equilibrate the constant pressure runs, 
while the constant volume simulations needed just 5,000 MCS for equilibration. 



cither close-packed structure as a reference crystal, which affords us a little more flexibility when it comes to choosing 
system sizes. The system sizes used are shown in Table [I] 

Once the crystal lattices have been chosen, one must choose a site mapping between crystals A and B . In the 
previous lattice-switch work for FCC and HCP structures more efficient switches which preserved close-packed planes 
were used. However, for binary-to-single-lattice switches, there are no such obvious similarities, therefore, we mapped 
sites between crystals at random. 

To predict the equilibrium phase behaviour, we must first determine the Gibbs free energy of all the candidate 
phases, as a function of the pressure, P* = Pa\/kT, and the composition, X = Nb/(Na + Nb), where Na and Nb 
are the numbers of A and B particles respectively. Naturally, we would seek to measure the Gibbs free energy directly 
via a lattice-switch simulation carried out in the constant pressure ensemble. However, when attempting to switch 
between the pure and the binary crystals using a single simulation cell, the fact that the equilibrium cell size is very 
different for these two structures means that the switch becomes ineffective. The lattice-switch only updates the lattice 
and the diameters, not the cell parameters, and therefore performing the switch from the reference lattice involves 
attempting to transform an equilibrium density pure crystal into a loosely packed binary crystal (as reducing the 
diameter of so many particles vastly increases the available volume). The conjugate move requires a transformation 
from an equilibrium density binary crystal to a pure crystal at such a high density that it is likely that all particles 
overlap. To facilitate either switching move, the volume of the simulation cell would have to be dilated so far that the 
system would rapidly melt, and indeed such a switch cannot be made to work in practice. It is possible, at least in 
principle, to overcome this limitation by adding an explicit volume dilation to the switch, so that the cell dimensions, 
the lattice and the diameters are all changed simultaneously. However, adding a volume transformation factor affects 
the probability of sampling the two phases in a non-trivial way. For example, a dilation to a large cell will almost 
always be accepted, whereas the compression to a small cell will be difficult to achieve. This will bias the measured 
free energy difference by a factor which could be determined by careful application of the mapping transformation 
matrix approach outlined in [8], but this introduces a number of computational overheads, and the slow switching 
rate mitigates against gathering good statistics. 

All these issues can be avoided in the constant volume ensemble by keeping the density fixed when switching between 
the reference and target structures, thus ensuring that the displacements are reasonable in both phases because the 
distances between particles are comparable in both phases under these conditions. 



B. Simulation Methodology 

To estimate the Gibbs free energy, three simulations were carried out for each state point. The first was a single 
phase constant pressure (NPT) simulation, used to estimate the equilibrium density and cell parameters (c/a ratio) 
for each binary crystal, at each pressure. Then, two lattice-switch NVT simulations were performed at the volume 
fractions and c/a ratios determined from the NPT simulations. The lattice switch attempts to transform between the 
target crystal (at the target density and c/a ratio) and a pure FCC crystal (at the same density, but with c/a = 1.0). 
The first run is used to evaluate the multicanonical weight function and the second uses this function to determine 



the free energy difference between the crystals (see section III C for more details) 
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Having determined the cell proportions, the density, and the Helmholtz free energy difference between the two 
structures (Ai* 1 ) for a range of pressures, we can then build the absolute free energy of the binary crystal by adding 
in the free energy of the pure phase. To determine this, we integrate the Alder equation of state [101 l4"Tl 142) . 



PV 3 

—— = + 2.566 + 0.55(V* - 1) - 1.19(F* - l) 2 + 5.95(V* - l) 3 , (10) 

NkT V* - 1 ; 

with respect to the reduced volume V* — V/Vq (where Vq is the close packed volume, Na 3 /^/2). This yields the 
Helmholtz free-energy: 



T^Alder / X/* i \ in r\A r 

-g— = -3 In — + 5.1241nF* - 20.78^* + -^-V* 2 - ^V* 3 + C 4 + 1, (11) 

NkT \ V* J 2 3 

where, following 42J, the constant of integration C4 is set to 15.05. The Gibbs free energy can be constructed from 
the Helmholtz form via: 



F 



Alder 



AF + PV - 1 - In 



PV 



NkT 



(12) 



where the last two terms arise from the difference between the Gibbs and Helmholtz free energies for an ideal gas. 
To compute the phase diagrams based on these free energies, we follow the approach of Bartlett [TH] and assume 
that the particles are immiscible in the solid phases, and that the binary fluid behaviour can be described accurately 
by the Mansoori equation [43 . At each pressure, the densities of the candidate phases are determined from the 
equations of state. From these densities, the Gibbs free energies and chemical potentials of the competing phases can 
be determined for each constant pressure line. The condition of equal chemical potentials is then applied along the 
constant pressure line, where the the coexistence between the crystal phases and the Mansoori fluid is estimated from 
the chemical potentials of each species in the fluid via (see [2U]): 



H f A + nfi f B = (l + n)xG(AB n ), (13) 

where \i A and /x^are the chemical potentials determined from the Mansoori fluid for the A and B particles respec- 
tively, and G(AB n ) is the Gibbs free energy of the solid phase that has n B particles for each A particle. 



C. Multicanonical Biassed Monte Carlo procedure 



In general, a lattice switch move from a zero-overlap microstate in one phase is unlikely to map onto a zero-overlap 
microstate in the other. Therefore, the simulation must monitor not only which phase the is currently being explored, 
but also how close the simulation is to being able to perform the lattice switch. To find a suitable order parameter, we 
first let M({u},A) denote the number of overlapping pairs of particles associated with the displacements {u}, within 
the structure A. We can then define the overlap order parameter: 



M({u}) = M({u},B)-M({u},A), (14) 

where M. is necessarily > (< 0) for realizable configurations of the A (B) structure. To ensure that the simulation 
can reach the so-called gateway states (A4 — 0), we apply the multi-canonical Monte Carlo method [44] as a way of 
biasing the simulation in a controlled fashion. To this end, we augment the system energy function such that 



£({u},A) = E({u} 7 A) + V (M({u})), (15) 

where r)(AA),AA = 0,±1,±2, ... constitute a set of multicanonical weights 44J. Given that a suitable weight 
function can be determined (see below), the canonical distribution P(A4\N, V) can be estimated from the measured 
multicanonical distribution P(A4\N, V, {77}) with the identification 
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v (N v) Em> P(M\N,V) E M >oP(M\N,V,{ V })e^ 
l ' j Em<oP(M\N,V) E M < P (M\N,V,{r)})en(M)- 

Here, the exponential reweighting maps the multicanonical distribution onto the canonical, unfolding the bias 
associated with the weights. 

The multicanonical weights are determined from non-equilibrium calculations, which are fast and give a fair ap- 
proximation to the equilibrium situation. We conduct a series of separate Monte Carlo simulations, where the particle 
displacements (and therefore M.) are reset to zero at the start of each. As these simulations will each relax to the 
equilibrium M. values, this forces the overall run to explore the full range of M.. Throughout this sequence of simu- 
lation, the program monitors and accumulates the macrostate transition probability matrix, as described in |45j . and 
at the end of the whole run, an estimate of the weight function is determined from the transitions probabilities by 
a simple shooting method |46j . This weight function is then used for the second NVT run, during which the time 
spent in each phase is recorded (with the biasing taken into account), and automatically block-analysed to produce 
an estimate of the free-energy difference. The actual timings varied depending on the target system and the number 
of particles, and are shown in Table [TJ 

The simulations operated under the usual Metropolis algorithm 47J, using random numbers generated by the 
Mersenne Twister algorithm |1S], as implemented in RngPack, version 1.1a [49 . For each move, a sphere is chosen 
at random, and a trial displacement is generated by taking the displacement of that particle Ui and adding a random 
three-dimensional vector Aui drawn uniformly from a cube of fixed size. The update in the displacement is accepted 
according to the usual Metropolis prescription: 



p a ({u} ^ {u} + Aui) = min{l,exp[-{£{{u + AiZ i },A)-£{{u},A))]}. (17) 

Therefore, moves that would lead to overlaps in the current phase are always rejected, and moves that would 
alter the number of overlaps in the conjugate phase are accepted or rejected based on the multicanonical weights. 
The magnitude of Aui was automatically adjusted to yield an acceptance rate of about one third, as this has been 
shown previously to produce an acceptable rate of decorrelation in A4 [46]. When M. = 0, the lattice-switch move is 
attempted once per Monte Carlo Sweep (MCS) (i.e. once every N particle moves) on average, and is always accepted 
as there is no energy cost or multicanonical weight difference associated with this move. 

Due to the short-range nature of the hard sphere interaction, the calculation of E({u},A) does not need to extend 
over all iV 2 particle pairs. To make the simulation more efficient, each particle only checks for interactions with its 
neighbours, up to a maximum interaction distance of 1.5 unit cells (i.e. 1.5 times the A-A separation). There are two 
lists of neighbours for each particle: one for each phase. As mentioned in section [TTJ we require a criterion to determine 
whether structure is stable, and the simulation achieves this by periodically scanning the system for particles that 
have moved further than the near neighbour distance away from their site (after having taken any centre of mass 
diffusion into account). If this is observed to happen at any point during any simulation at a given pressure, then the 
structure is considered unstable, and no longer considered as a candidate structure at that pressure. 

To establish the equations of state for these systems, a number of constant pressure simulations were also performed. 
In this case, the simulation cell parameters (the lengths of the sides of the cell, and therefore the volume) can fluctuate 
during the simulation. Such moves are accepted with probability: 



Pa(V -» V) = min{l,exp[-(£ ({u}',A) - £({u},A)) - P*AV + Nbi(V'/V)]}, (18) 

where V' is the volume associated with the trial move. Of course, this cell-size dilation can change £ globally, and 
so an expensive energy calculation must be carried out for each trial volume move. Volume moves were attempted at 
random, once per MCS on average. The maximum size of the random volume-dilation is automatically adjusted so 
that the acceptance rate for volume moves is about one half, which produces an acceptable rate of exploration of the 
volume parameter space |4"6] . 



IV. RESULTS 



A. Pressure curves and stability 

The equations of state, as measured by Monte Carlo simulation in the constant pressure ensemble, are shown in 
Figure [2] The agreement with the data presented in [7] indicates that the Molecular Dynamics simulation presented 
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Figure 2: Equations of state determined in the constant pressure ensemble, measuring the density as a function of the pressure. 
Errors in the density estimates are smaller than the symbol size. The Alder equation of state for the pure crystal [40] is shown 
as a solid line, and the equation of state for AB13 as observed in [7] is shown as a solid gray line. 



there used the smaller (N — 112) system size when computing the equation of state. Also, the close agreement 
between the Alder equation of state [40 and the simulated FCC crystal is clear. However, we note that while the 
Alder form (eq. 10 1 provides a reliable estimate of the equation of state, the individual terms are difficult to reproduce. 
Despite having determined the density of the FCC structure at each pressure to within at least 0.1%, attempting to 
fit an equation of the Alder form to this data led to a poorly-determined fit that only reproduced the original terms 
to within one significant figure. It is therefore clear that a much larger range of pressures must be used to determine 
Alder-style equations of state for these structures, and given that quite a modest range of pressures was explored in 
the original work |41j . it is not clear that the parameterization supplied in cq. 10 is a particularly accurate one. 

During the equation of state simulations, the cell proportions were also monitored and Figure [3] shows the results 
in terms of the c/a ratio. As expected on symmetry grounds, the AB13 and CsCl crystals are stable at c/a = 1. 
However, the AB2 crystal tends to expand in the direction perpendicular to the alternating stacking planes of A and 
B particles, and the degree of expansion found here is consistent with the results presented in [2]. The c/a ratio 
decreases as the pressure is lowered towards the fluid region, and as the simulations approach the melting point, the 
fluctuations in the volume and aspect ratio of all the crystals become large as the structures start to break down. 
As noted earlier, this break-down condition discounts extended metastability and the possibility that the structure 
may be stabilised by some equilbrium concentration of defects, and so provides an upper bound on the actual limit 
of stability of the structure. 



B. Free energies 

The free energy difference between the reference (FCC or HCP) and target crystals AB2 and AB13 at a = 0.58, 
and for AB(CsCl) at a — 0.73 are shown in Figure [4] In all cases, the free energy difference between FCC and HCP 
1 x 10 _3 fcT per particle) is so small that it cannot be resolved on the scale of these plots. The AB2 and AB13 
structures are lower in free energy than the FCC crystals for all densities at these size ratios. 

Figure [4] also shows that the AB(CsCl) structure has a significantly higher Helmholtz free energy than FCC at high 
densities, the difference falls rapidly as the density is lowered towards melting. However, it is important to remember 
that the AB(CsCl) structure is not competing with a pure FCC crystal, but with one or more coexisting phases at 
equal pressures and chemical potentials. Figure [5] illustrates this by comparing the Gibbs free energy of AB(CsCl) 
with those of the competing phases - a Mansoori fluid, a pure FCC crystal coexisting with a B-rich Mansoori fluid, 
and a dual crystal state made up of separate A and B FCC crystals. The AB(CsCl) structure is clearly not stable 
over this range, and this result compares favourably with results from DFT calculations in [31] and via cell theory in 
[29j . However, we note that the free energy difference is quite small, and so a small amount of charge asymmetry may 
be sufficient to stabilise the structure in the low density solid phase. 

For all structures, the free energy difference measured in the constant density ensemble shows very little dependence 
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Figure 3: Observed c/a ratios for the target phases. 
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Figure 4: The raw Helmholtz free energy difference per particle between the reference FCC and target AB2 and AB13 crystals 
at q = 0.58, and AB(CsCl) at a = 0.73, as measured by the lattice switch simulations. The final data points at low <j> indicate 
the lowest density at which the structure was determined to be stable. The points at the high density end indicate the upper 
limit of the chosen pressure range explored by the simulations. Some simulations were performed at larger system sizes, as 
shown here, and indicated that the overall form of the free energy difference does not depend strongly on the system size. 



on system size, and even rather small systems provide reasonable estimates of the free energy difference with respect 
to FCC. To allow comparison with the literature results, the Helmholtz free energies of the structures for a = 0.58 
are shown in figure [6j and the Gibbs free energies (computed via eq. |12[ ) are shown in Figure [7] The lattice switch 
results are in clear agreement with the results determined via integration methods [6j [7] , and with the DFT results 
presented in J3UJ and the cell theory results in [25] . For clarity, Figure [6] only includes the results presented in [5] [7] 
and omits the almost identical results of [25] [25] . 
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Figure 5: Gibbs free energies as a function of pressure for AB(CsCl) a = 0.72 at a fixed composition of X = 0.5. The Gibbs 
free energies of the Mansoori fluid, the pure A FCC crystal coexisting with a Mansoori fluid, and the dual FCC A and B crystal 
phases are also shown. Finite size effects cannot be resolved on this plot. 




Figure 6: The excess (i.e. with respect to an ideal gas) Helmholtz free energies of AB2, AB13 and FCC at a — 0.58. Data from 
literature for AB2 |B] and AB13 [7] is also shown for comparison. 



C. Phase diagrams 



Following the prescription outlined in section |IIIB| the lattice switch data were used to predict the pressure- 
composition phase diagram for the binary hard sphere system at a = 0.58, and the results are shown in Figure [8] 
The corresponding partial volume-fraction phase diagram is shown in Figure [9j Both these figures agree well with the 
established results [SJ US] , although the partial volume fraction phase diagram does not extend to as high densities as 
the one published in [5 due to the limited range of pressures explored in the current work. 
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Figure 7: Gibbs free energies for AB2 and AB13 at a — 0.58. The Gibbs free energies of Mansoori fluid and the A FCC crystal 
coexisting with a B-rich Mansoori fluid at the appropriate compositions are also shown. Finite size effects cannot be resolved 
on this plot. 
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Figure 8: P-X Phase Diagram at a = 0.58, determined from the lattice switch results using the methodology described in 
section 1111 Bl 



V. DISCUSSION 



We have presented a method by which the lattice switch approach can be applied to binary mixtures, and applied 
it to binary hard sphere system. The technical innovation introduced here involves including a rescaling of particle 
sizes in the switch. Combining our free energy differences enables us to build a phase diagram which is in excellent 
agreement with previous theoretical results at the size ratio 0.58. The only errors in the lattice-switch MC approach 
come from finite size effects and sampling errors, which are easy to quantify. The phase diagrams themselves rely on 
the Alder and Mansoori equations of state for the FCC and liquid free energies respectively. 

We also investigated the experimentally reported CsCl phase at size ratio 0.73, and found it to be unstable compared 
to phase-segregated FCC or the liquid. However the energy differences are small, and it is possible that either charge 
effects or polydispersity may be responsible for the experimental observation of this structure. 

The dependence upon the Alder crystal is perhaps unfortunate, as it is unclear how accurate the parameters of this 
equation of state are. However, the Alder equation of state is dominated by the leading term and rather insensitive 
the poorly-determined parameters. Our results are in good agreement with the literature, including those [5] which do 
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not depend on the Alder formulation (due to using an Einstein crystal as a reference state). In the same vein, building 
a lattice-switch mapping to an analytically soluble reference crystal would also allow the use of the Alder equation of 
state to be avoided. This Hamiltonian-switch approach could map from an Einstein crystal or single-occupancy cell 
model to the target potential and structure, and would allow absolute free energies to be measured directly. 

Other future work includes lower size ratio, although as mentioned before, the mobile smaller spheres will present 
challenges, although it may be possible to use tethers, in a manner similar to that employed for the fluid-solid 
transition [10]. Also, these structures have been observed for noble gas mixtures [23 , and this could be investigated 
by generalising this approach in the same manner as the extension to soft-potentials of the original lattice switch 
approach [TT] . 
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